#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

   //////////////////
   // Input   ///////
   //////////////////
   char dir[1024];
   int Nmass;
   int  set_p[4];
   char setpprint[16];
   char Oa2print[16];
   int iprop;
   int beta100;

   char momentum_order[1024];
   ////// Parsing command line

   for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
	    if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
	    if(strncmp(argv[i],"-beta100",8)==0)  sscanf(argv[i],"-beta100=%d",&beta100); 
            if(strncmp(argv[i],"-set_p",6)==0)  sscanf(argv[i],"-set_p=%d,%d,%d,%d",&set_p[3], &set_p[2], &set_p[1], &set_p[0]);
	
	    if(strncmp(argv[i],"-iprop",6)==0)  sscanf(argv[i],"-iprop=%d",&iprop);
	    if(strncmp(argv[i],"-Oa2",4)==0)  sscanf(argv[i],"-Oa2=%s",Oa2print);

             if(strncmp(argv[i],"-momentum_order",15)==0)  sscanf(argv[i],"-momentum_order=%s",momentum_order);

  }

   int Nmass_comb=Nmass;
   for (int i=1; i<Nmass; i++)  Nmass_comb=Nmass_comb+(Nmass-i);
   /////////////////////////
  //// Momenta set
  //////////////////////////////

  whichProp(iprop,setpprint);
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom_%d.dat",beta100);
  read_momenta_data(momentafile);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(momentum_order,nn);
  Momentum_Magro(momentum_order,set_p);
  // equivalent momenta
  Nmom_eq=Number_of_Equivalent_Momenta(SinP2);
  Average_over_Equivalent_Momenta(SinP2,SinP2,SinP2_eq);

  //////////////////////////////////
  // IO files
  ///////////////////////////////

  char file[1024];
  char dir_qav[1024];
   char pm[8];

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	//////////////////////////////////
	// Read Zs
       // NOTE: Lambda & Zq MUST be reals (check in the last step)
       // I work with reals from now on
	/////////////////////////////



	double *****ZO=malloc(sizeof(double****)*(3));
        for (int ipm=0; ipm<3; ipm++) ZO[ipm]=malloc(sizeof(double***)*(5));
	for (int ipm=0; ipm<3; ipm++) for (int igamma=0; igamma<5 ; igamma++)  ZO[ipm][igamma]=malloc(sizeof(double**)*(Nmass_comb));
        for (int ipm=0; ipm<3; ipm++) for (int igamma=0; igamma<5 ; igamma++)  for (int imass=0; imass<Nmass_comb; imass++) ZO[ipm][igamma][imass]=malloc(sizeof(double*)*(Nmom_eq));
        for (int ipm=0; ipm<3; ipm++) for (int igamma=0; igamma<5 ; igamma++) for (int imass=0; imass<Nmass_comb; imass++) for (int imom=0; imom<Nmom_eq; imom++) ZO[ipm][igamma][imass][imom]=malloc(sizeof(double)*(Nclust+1));

        double ****Zq=malloc(sizeof(double***)*(3));
	  for (int ipm=0; ipm<3; ipm++) for (int imass=0; imass<Nmass; imass++) Zq[ipm]=malloc(sizeof(double**)*(Nmass));
        for (int ipm=0; ipm<3; ipm++) for (int imass=0; imass<Nmass; imass++)  Zq[ipm][imass]=malloc(sizeof(double*)*(Nmom_eq));
        for (int ipm=0; ipm<3; ipm++) for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++) Zq[ipm][imass][imom]=malloc(sizeof(double)*(Nclust+1));

        double *aux_clust=malloc(sizeof(double)*(Nclust+1));

    // Zq
    for (int imass=0; imass<Nmass; imass++) {
    for (int ipm=0; ipm<2; ipm++)
    {
    whichpm(ipm,pm);
    sprintf(file,"Zq_%s_m%d",Oa2print,imass);
    sprintf(dir_qav,"%s%s/%s/QAV/",dir,pm,setpprint);
    for (int imom=0; imom<Nmom_eq; imom++)
    {
		read_input_onemomenta(dir_qav,file,aux_clust,imom);
    		for (int iclust=0; iclust<Nclust+1; iclust++) Zq[ipm][imass][imom][iclust]=aux_clust[iclust];
    }
    }

    }
   // prom
   for (int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) Zq[2][imass][imom][iclust]=(Zq[0][imass][imom][iclust]+Zq[1][imass][imom][iclust])/2.;
    

   char gammaprint[8];
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
        int imass=0; 
        for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
        {
	sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
	    for (int ipm=0; ipm<2; ipm++)
	    {
	    whichpm(ipm,pm);
	    sprintf(dir_qav,"%s%s/%s/QAV/",dir,pm,setpprint);
	    for (int imom=0; imom<Nmom_eq; imom++)
	    {
                read_input_onemomenta(dir_qav,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) ZO[ipm][igamma][imass][imom][iclust]=aux_clust[iclust];
            }
	    }
	imass++;
	}
    }
   for (int imass=0; imass<Nmass_comb; imass++) for (int igamma=0; igamma<5; igamma++) for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) ZO[2][igamma][imass][imom][iclust]=(ZO[0][igamma][imass][imom][iclust]+ZO[1][igamma][imass][imom][iclust])/2.;


   // print output
   sprintf(dir_qav,"%sa/%s/QAV",dir,setpprint);
   for (int imass=0; imass<Nmass; imass++) {

         sprintf(file,"Zq_%s_m%d",Oa2print,imass);
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]= Zq[2][imass][imom][iclust];
                 write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq,imom);
          }
          sprintf(file,"%s/Z_ALLP/Zq_%s.dat",dir_qav,Oa2print);
          FILE *fout=fopen(file,"w");
          if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
          }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Zq[2][imass][imom][iclust];
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],Zq[2][imass][imom][Nclust],error_jack_real(aux_clust,Nclust));
        }
      fclose(fout);

   }
   for (int igamma=0; igamma<5; igamma++)
   {
        whichGammaProm(igamma,gammaprint);
	int imass=0;
	for (int im2=0; im2<Nmass; im2++) for (int im1=im2; im1<Nmass; im1++)
	{
        sprintf(file,"Z%s_%s_m%d%d",gammaprint,Oa2print,im1,im2);
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]= ZO[2][igamma][imass][imom][iclust];
                 write_output_onemomenta_real(dir_qav,file,aux_clust,SinP2_eq,imom);
        }
          sprintf(file,"%s/Z_ALLP/Z%s_%s.dat",dir_qav,gammaprint,Oa2print);
          FILE *fout=fopen(file,"w");
          if(fout == (FILE *)NULL) {
              fprintf(stderr,"Unable to open file %s for writing\n",file);
              return EXIT_FAILURE;
          }
         for (int imom=0; imom<Nmom_eq; imom++) {
                for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=ZO[2][igamma][imass][imom][iclust];
                fprintf(fout,"%8.10f %8.10f %8.10f  \n",SinP2_eq[imom],ZO[2][igamma][imass][imom][Nclust],error_jack_real(aux_clust,Nclust));
        }
      fclose(fout);
	imass++;
	}

   }

   return 0;

}
